# https://zhuanlan.zhihu.com/p/367213524


import cv2
import matplotlib.pyplot as plt
import numpy as np
from skimage.feature import graycomatrix, graycoprops
import copy

import detector

def calc_glcm_coprops(img, graylevel,dist,degree,show=False,main_title=None):
    # 压缩图片灰度级
    tablelevel = np.array([(i//graylevel) for i in range(256)]).astype("uint8")  
    img_graylevel = cv2.LUT(img, tablelevel)  # 灰度级压缩为 [0,graylevel-1]
    glcm = graycomatrix(img_graylevel, dist, degree, levels=graylevel)  
    
    # 这个字典按照统计量存放了所有不同距离与角度下计算出来的灰度共生矩阵的统计量
    # 每个键所对应的值都是一个形状为 [len(dist), len(degree)] 的矩阵
    coprops_by_Statistics = {'contrast':0,
                            'dissimilarity':0,
                            'homogeneity':0,
                            'energy':0,
                            'correlation':0,
                            'ASM':0
    }
    # 逐个统计量解出
    for prop in coprops_by_Statistics.keys():
        coprops_by_Statistics[prop] = graycoprops(glcm, prop).round(8)

    # 这个字典作为元素模板，用于创建字典数组，且用于暂时存放单个 glcm 的所有统计量
    coprops = {'contrast':0,
                'dissimilarity':0,
                'homogeneity':0,
                'energy':0,
                'correlation':0,
                'ASM':0
    }
    # 先创建数组的行，行中的字典属于相同距离，不同角度的glcm
    # 这个一维字典数组作为二维字典数组的元素
    coprops_by_glcm_same_dist = []
    for j in range(len(degree)):
        coprops_by_glcm_same_dist.append(copy.deepcopy(coprops))

    # 这个字典按照不同灰度共生矩阵存放了所有不同距离与角度下计算出来的灰度共生矩阵的统计量
    # 形状为[len(dist), len(degree)]，每个元素都是一个字典
    coprops_by_glcm = []
    for i in range(len(dist)):
        coprops_by_glcm.append(copy.deepcopy(coprops_by_glcm_same_dist))

    for i in range(len(dist)):
        for j in range(len(degree)):
            for prop in coprops.keys():
                coprops[prop] = coprops_by_Statistics[prop][i,j]
            coprops_by_glcm[i][j] = copy.deepcopy(coprops)

    if show:
        plt.figure(figsize=(len(degree)*3, len(dist)*3))
        plt.suptitle(main_title)
        for i in range(len(dist)):
            for j in range(len(degree)):
                plt.subplot(len(dist),len(degree),i*len(degree)+j+1)
                plt.title(r"d={},$\theta$={:.2f}".format(dist[i], degree[j]))
                plt.text(1,12,get_coprops_fstring(coprops_by_glcm[i][j]),color='green',size=14)
                plt.imshow(glcm[:,:,i,j], 'gray')
        plt.tight_layout()
        detector.display_imgs([img_graylevel], ["img_graylevel"],num_cols=1,main_title=main_title)
    return glcm, coprops_by_glcm

def get_coprops_fstring(coprops):
    coprops_fstring = "contrast:{}\n".format(coprops['contrast'])
    coprops_fstring += "dissimilarity:{}\n".format(coprops['dissimilarity'])
    coprops_fstring += "homogeneity:{}\n".format(coprops['homogeneity'])
    coprops_fstring += "energy:{}\n".format(coprops['energy'])
    coprops_fstring += "correlation:{}\n".format(coprops['correlation'])
    coprops_fstring += "ASM:{}".format(coprops['ASM'])
    return coprops_fstring


# 14.11 特征描述之灰度共生矩阵 (skimage)
if __name__ == "__main__":
    single_img_path = "../fabric-defect/油污+缺纬/T03903_00.bmp"
    img = detector.cv_read(single_img_path)
    img_EquHisto = detector.AdaptiveEquHisto(img)

    # 取左上角的ROI
    top_left = [0,0]
    buttom_right = [120,200]
    roi_top_left = detector.get_roi(img_EquHisto,top_left, buttom_right)

    # 取右上角的ROI
    [height, width] = img_EquHisto.shape

    top_left = [width-121, 0]
    buttom_right = [width-1,200]
    roi_top_right = detector.get_roi(img_EquHisto,top_left, buttom_right)

    ksize = 100 
    params = {'ksize':(ksize, ksize), 
                    'sigma':5, 
                    'theta':0,                                # theta代表条纹旋转角度
                    'lambd':8,                                   # lambd为波长 波长越大 条纹越大
                    'gamma':1,                                  # gamma越大核函数图像越小，条纹数不变，sigma越大 条纹和图像都越大
                    'psi':0,                                   # psi 是相位偏移 psi这里接近0度以白条纹为中心，180度时以黑条纹为中心
                    'ktype':cv2.CV_32F}

    kern = cv2.getGaborKernel(**params)
    roi_top_left_gabor = cv2.filter2D(roi_top_left, -1, kern)
    roi_top_right_gabor = cv2.filter2D(roi_top_right, -1, kern)

    slope_top_left = detector.calc_mean_slope(roi_top_left_gabor)
    slope_top_right = detector.calc_mean_slope(roi_top_right_gabor)

    '''
    已知斜率k与一点(x0, y0)，又知另一点的一个坐标y, 求直线方程上另一点的x。

    (y - y0) = k*(x - x0)
    y - y0 = k*(x - x0)
    (y - y0)/k = x - x0
    x = (y - y0)/k + x0
    '''
    # 透视变换，需要取原图的4个顶点，变换到新的4个位置上
    point_top_left = [0,0]
    point_top_right = [width-1, 0]
    point_buttom_left = [0, height]
    point_buttom_right = [0, height]

    point_buttom_left[0] = (point_buttom_left[1] - point_top_left[1]) / slope_top_left + point_top_left[0]
    point_buttom_right[0] = (point_buttom_right[1] - point_top_right[1]) / slope_top_right + point_top_right[0]

    points = [point_top_left, point_top_right, point_buttom_left, point_buttom_right]
    img_draw_point = detector.draw_points(img_EquHisto,points,(255,0,0),3,8)

    dst_width = point_buttom_right[0]-point_buttom_left[0]
    dst_width = np.round(dst_width).astype(int)
    dst = [[0,0],[dst_width,0],[0,height],[dst_width,height]]
    img_Standardize = detector.PerspectiveTransform(img_EquHisto,points,dst,(dst_width,height))
    detector.cv_show("img_Standardize", img_Standardize)

    # 压缩图片灰度级
    # height, width = img_Standardize.shape
    table16 = np.array([(i//16) for i in range(256)]).astype("uint8")  # 16 levels
    gray16 = cv2.LUT(img_Standardize, table16)  # 灰度级压缩为 [0,15]

    # 计算灰度共生矩阵 GLCM
    dist = [1, 70, 140, 210]  # 计算 4 个距离偏移量 
    degree = [0, np.pi/4, np.pi/2, np.pi*3/4]  # 计算 4 个方向
    glcm = graycomatrix(gray16, dist, degree, levels=16)  # 灰度级 L=16
    print(glcm.shape)  # (16,16,4,4)

    # 由灰度共生矩阵 GLCM 计算特征统计量
    for prop in ['contrast', 'dissimilarity','homogeneity', 'energy', 'correlation', 'ASM']:
        feature = graycoprops(glcm, prop).round(4)  # (4,4)
        print("{}: {}".format(prop, feature))

    plt.figure(figsize=(9, 6))
    plt.suptitle("GLCM by skimage, youcans")
    for i in range(len(dist)):
        for j in range(len(degree)):
            plt.subplot(len(dist),len(degree),i*len(degree)+j+1)#, plt.axis('off')
            plt.title(r"d={},$\theta$={:.2f}".format(dist[i], degree[j]))
            plt.imshow(glcm[:,:,i,j], 'gray')
    plt.tight_layout()
    plt.show()
# ————————————————
# 版权声明：本文为CSDN博主「youcans_」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
# 原文链接：https://blog.csdn.net/youcans/article/details/125693533